#ifndef GSL_MATRIX_CACULATE_H
#define GSL_MATRIX_CACULATE_H

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_linalg.h>

class  GSLMatrixCaculate
{
public:

	GSLMatrixCaculate();
	~GSLMatrixCaculate();
	void GetInverse(gsl_matrix *A, gsl_matrix *inverse);
	gsl_matrix MatrixMultiply(gsl_matrix *A, gsl_matrix *B, gsl_matrix *C);
	double GSL_Vector_Mean(gsl_vector *A);

private:

};

GSLMatrixCaculate::GSLMatrixCaculate()
{

}
GSLMatrixCaculate::~GSLMatrixCaculate()
{

}
void GSLMatrixCaculate::GetInverse(gsl_matrix *A, gsl_matrix *inverse)
{
	int n = A->size1;
	gsl_matrix *tmpA = gsl_matrix_alloc(n, n);
	gsl_matrix_memcpy(tmpA, A);
	gsl_permutation *p = gsl_permutation_alloc(n);
	int sign = 0;
	gsl_linalg_LU_decomp(tmpA, p, &sign);
	//inverse = gsl_matrix_alloc(n, n);
	gsl_linalg_LU_invert(tmpA, p, inverse);
	gsl_permutation_free(p);
	gsl_matrix_free(tmpA);
}

gsl_matrix GSLMatrixCaculate::MatrixMultiply(gsl_matrix *A, gsl_matrix *B, gsl_matrix *C)
{
	size_t rowA = A->size1;
	size_t colA = A->size2;
	size_t rowB = B->size1;
	size_t colB = B->size2;
	if (colA != rowB)
	{
		//GSL_ERROR("matrix must have same Col and Row", GSL_EBADLEN);
	}

	else
	{
		//if (c.data) delete[] c.data;
		//C = gsl_matrix_alloc(rowA,colB);
		//c.data = new double[c.rows * c.cols];
		for (size_t i = 0; i<rowA; i++)
		{
			for (size_t j = 0; j<colB; j++)
			{
				double sum = 0;
				for (size_t k = 0; k<colA; k++)
				{
					sum = sum + (gsl_matrix_get(A, i, k)) * (gsl_matrix_get(B, k, j));
				}
				gsl_matrix_set(C, i, j, sum);
				//C->data[i, j] = sum;
				//cout << C->data[i, j]<<endl;
			}
		}
	}
	return *C;
}

double GSLMatrixCaculate::GSL_Vector_Mean(gsl_vector *A)
{
	int sizeA;
	double VectorMean;
	double sum = 0;
	sizeA = A->size;
	for (int i = 0; i < sizeA; i++)
	{
		sum += gsl_vector_get(A, i);
	}
	VectorMean = sum / sizeA;
	return VectorMean;
}
#endif